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We theoretically study the stability of lattice supersolid states in the extended Bose-Hubbard 
model with bounded spatial disorder. We construct a disorder mean field theory and compare with 
quantum Monte Carlo calculations. The supersolid survives weak disorder on the simple cubic 
lattice. We also find that increasing disorder strength can transform a lattice solid into a supersolid 
as it tends to percolate through the disorder landscape. 
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Studies of the extended Bose-Hubbard (BH) model 
seek to capture the essential properties of a wide variety 
of physical systems including helium [TJ [2J, Josephson 
junction arrays [3. , and certain narrow-band supercon- 
ductors [4: . The model reveals a lattice supersolid phase 
(a simultaneous superfluid and solid) in its mean field 
phase diagram [2j [3]. Recent experiments in solid 4 He 
demonstrate evidence [5] for the long-anticipated super- 
solid phase. Disorder and defects are believed to impact 
the robustness of the observed phase [B] . 

Theoretical studies of the lattice supersolid suggest 
that it is rather delicate. For example, quantum Monte 
Carlo (QMC) analyses reveal that large interaction 
strengths are required to stabilize the supersolid against 
quantum and thermal fluctuations in low dimensions [3- 
[5J . More troubling results on the square lattice [TO] show 
that a different type of perturbation, external spatial dis- 
order, destroy the solid itself leaving no chance for the 
supersolid. The latter results stem from an Imry-Ma- 
type argument [101 E] implying that the solid is unsta- 
ble in the presence of arbitrarily weak disorder in less 
than three dimensions (3D). Recent QMC results firmly 
established that in 3D, e.g., on the simple cubic (SC) 
lattice [12] . the supersolid is in fact stable against both 
thermal and quantum fluctuations yet, to the best of our 
knowledge, no QMC studies on the effect of spatial disor- 
der on the SC lattice supersolid have been reported. On 
the SC lattice we expect the stability of the solid against 
disorder [5] but the fate of the supersolid is unknown. 

Disorder in BH models generates an intriguing com- 
pressible glassy state, the Bose-glass [13], that could 
compete with supersolids. Considerable controversy (see 
Ref. [14] and references therein) regarding the precise lo- 
cation of the Bose-glass in phase diagrams of the ordi- 
nary BH model (the short-range limit of the extended 
BH model) with disorder stems from Griffiths effects. 
Griffiths effects arise from statistically rare, but relevant, 
disorder configurations that can generate glassy states. 
These rare configurations combine with strong quantum 
fluctuations over large length scales. As a result they 
are often missed in mean field theories (MFTs) and can 
mislead finite-sized QMC studies. 

The theorem of inclusions [T3] establishes a key feature 



in the phase diagram topology of disordered BH models: 
We expect a compressible glassy state to completely sur- 
round incompressible states. But in the extended BH 
model, supersolids are also expected to lie near the in- 
compressible solid in parameter space. The theorem of 
inclusions thus mandates a competition between glassy 
states and supersolids in extended BH models, opening 
the possibility that supersolids do not survive even weak 
disorder on cubic lattices. 

Work on the relationship between quantum order and 
disorder has become more pressing with the advent of 
cold atomic gas experiments demonstrating controlled 
disorder [T5HT8] . Disorder in optical lattice experi- 
ments, imposed with speckle laser light or bichromatic 
incommensuration, exhibit intriguing insulating phases 
of bosons [TMPE] . Long range interactions among bosons, 
caused by dipolar moments [T5], band effects [20], or 
other mechanisms [2T] suggest that the disordered ex- 
tended BH model could underlie the essential properties 
of future optical lattice experiments. 

We examine the stability of disordered lattice super- 
solids in the presence of quantum and thermal fluctua- 
tions. We construct a MFT phase diagram of the soft- 
core extended BH model with disorder and compare with 
QMC. For weak disorder we find that, in contrast to hard- 
core bosons on the square lattice [TUJ, the solid and the 
supersolid are indeed stable on the SC lattice. We then 
increase disorder to examine the role of strong disorder 
on bosonic solids. We find a striking effect: Disorder 
transforms the solid into a percolating supersolid in di- 
rect analogy to disorder triggered superfluidity via per- 
colation in BH models [14] [22H24] . We construct a mean 
field site-percolation picture to qualitatively capture the 
disorder triggered supersolid behavior. 

We study the interplay between disorder and super- 
solids with the extended BH model: 

H = 1 E { b h + H - c -) + f E n i( n J x ) 

(id) 3 

where the operator 6] creates a boson at the lattice site 
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Table I. List of possible phases and related order parameters: 
static structure factor SV, superfluid stiffness p s , and com- 
pressibility K. 
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j and the number operator is rij = bjbj. The first term 
imposes a hopping energy gain (t) among nearest neigh- 
bors. The second and third terms impose energy penal- 
ties for multiple occupancy at a site (U) and occupancy 
of nearest neighbor sites (V), respectively. We consider 
a bipartite lattice with periodic boundaries. To study 
a regime consistent with spatially decaying interactions 
and a strong supersolid we choose zV = U = 1 |5J [12], 
where z is the coordination number. Our MFT applies to 
any bipartite lattice but when using QMC with Eq. ^ 
we will work on the SC lattice, z = 6, with L 3 sites. 
The last term denotes a spatially varying chemical po- 
tential: [ij = fi + 5[ij with bounded disorder defined by 
the random number: b\ij G [—A, A]. 

We first examine the low-temperature phase diagram 
with a MFT that includes spatial disorder and thermal 
fluctuations but excludes quantum fluctuations. We de- 
couple sites by defining two density and two phase or- 
der parameters: m a = (n a ) and <j) a = (b^), respec- 
tively a € {A } B} denotes the sublattice index. Here 
and in the following the expectation values include disor- 
der averaging. The mean field Hamiltonian |12j becomes: 
H MF = Ha+Hb + C, where: 

H a = -zt{tf a + b a )(f> a > + ^n a (n a - 1) 

+ zVn a m a i ~ p,n a , (2) 

and C = 2zt4>A<t>B — zVuiaitib- In the above the prime 
denotes a ^ a' . fi = fi+e contains a continuously tunable 
parametrization of disorder, e. 

Self-consistent minimization of the free energy with re- 
spect to the order parameters solves the mean field equa- 
tions. The mean field free energy is: 

F MF = C - P~ 1 In Tr { e -^+H B ) j 

= C-r 1 \nJ2e-^, (3) 

a, 7 

where /? is the inverse temperature and E an denotes the 
7 th energy eigenvalues of H a . In the following we work 
at low densities. We find that a restricted Fock num- 
ber eigenstate basis accurately captures the low density 
properties: \N) — {|0), |1), |2)}. We then diagonalize the 



3x3 matrix (N\H a \N') to obtain three eigenvalues for 
each a, E an , where 7 = 1,2, 3. 

We use the free energy to solve the coupled mean field 
equations. We numerically solve the following four in- 
tegral equations: 7[<9^F MF ] = and I[d„ la F MF ] = 0, 
for m a and </> Q . Here I[R] = J^ A de(R/2A) denotes an 
integral over bounded disorder. 

The mean field equations yield several different ground 
states that are defined by combinations of disorder- 
averaged order parameters. Solid order is defined by 
long-range oscillations in the density-density correlation 
function (diagonal long-range order in the density ma- 
trix) or, equivalently, peaks in the static structure factor 
at wavevector, Q: 

S Q ^L- 6 5>^- rfc >(n,n fc ), (4) 

that indicate a spontaneous breaking of the sublattice 
symmetry. For the large values of V considered here an 
oscillation of the density between sublattices is favored, 
Q = (7r,7r,7r) on the SC lattice. In the mean field limit 
this corresponds to S^ F = (rriA — m B ) 2 /A. 

The superfluid phase is defined by off-diagonal long 
range order in the density matrix. The superfluid stiff- 
ness is: 

p^L- 6 ^(6]6 fe + H.c). (5) 

3,k 

This representation is well approximated by the stiffness 
measured in QMC, (W 2 )/3t/3, where W is the winding 
number [T2J [55] . In the mean field limit Eq. ^ yields 
^mf = ^ A _|_ B )2y2. The supersolid is defined by co- 
existing superfluid and solid order. The first three rows 
of Table [I] indicate expected orders in the uniform limit, 
A = 0. 

New phases arise in the presence of disorder, A > 0. 
For V = U/z we expect the solid phase to dominate 
but at low [i the spatially varying part of the chemi- 
cal potential can disrupt the perfect solid. Here even 
weak disorder can create compressible domains that ex- 
tend across the system. We define a disordered solid as 
a solid which has a finite S^, but zero global compress- 
ibility, k = (3((n 2 ) - (n) 2 ) (Table 

We first use MFT to test the stability of the supersolid 
against weak disorder. The MFT phase diagram (Fig. [lj 
shows that the solid and the supersolid remain stable 
for intermediate chemical potentials A < /i < Vz — A. 
(The top panel of Fig. [2] shows MFT order parameters 
for a characteristic line in the phase diagram at non-zero 
temperature.) But for /1 < A (or n > Vz — A), disorder 
will tend to push the chemical potential beyond the gap 
of the uniform solid. The energy gap (at A = 0) for the 
solid is given by the width of the solid lobe at fixed t, 
Eg — Mmax — Mmini m the /i vs. t phase diagram. Here 
Mmax (Mmin) i s the largest (smallest) chemical potential 
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Figure 1. (Color online) Solid lines show the zero temperature 
phase diagram of the extended Bose-Hubbard model, Eq. |l]), 
with weak disorder, A/U = 0.1, computed with MFT. The 
strong nearest neighbor repulsion zV/U = 1 opens a large 
solid lobe surrounded by the supersolid. The phase diagram 
is similar to the non disordered case except for low fj, (below 
the horizontal solid line) where the compressibility is finite. 
The dotted line marks a possible transition between different 
supersolids. The dashed line is a schematic of the location of 
an expected quantum Griffiths- type region. 



allowing the solid phase at fixed t. 

Fig. [I] shows that for fi < A the solid becomes com- 
pressible, k > 0. The horizontal line located at /x = A 
marks a transition from the solid to the disordered solid. 
(At finite temperatures this line becomes a crossover in 
MFT.) The MFT is defined locally and therefore just 
measures a local compressibility. In what follows we will 
describe QMC calculations. Our preliminary QMC re- 
sults show k > for fj, < A. We tentatively conclude 
that the disordered solid appears for it < A. 

The dotted line separates two different supersolid re- 
gions. MFT shows a cusp in the stiffness along the dotted 
line. Below the dotted line we expect a supersolid defined 
by combining a superfluid with a disordered solid that 
has compressible domains with a different spatial scal- 
ing behavior than above the dotted line. Domain size 
scaling and transitions marked by horizontal lines will be 
explored in future work. In the grand canonical ensemble 
phase separation (between a superfluid and a solid) man- 
ifests as a density discontinuity along phase boundaries 
[51 |H] . Our preliminary QMC results indicate phase sep- 
aration only near the dashed line for t/U ~ 0.04 — 0.05. 

We study the impact of quantum fluctuations using 
QMC. Equation ([I]) does not have a sign problem and 
therefore allows numerically exact studies using QMC 
on finite sized systems. We use the stochastic series ex- 
pansion representation with directed loop updates [26] 
within the Algorithms and Libraries for Physics Simula- 
tions (ALPS) framework [37J to evaluate order parame- 
ters of Eq. ([!]). We average over ~ 500 — 1000 disorder 
realizations for each data point reported. For each dis- 
order configuration a sufficient number of thermalized 
QMC steps are chosen to yield convergence of QMC er- 
ror. Error bars indicate the standard deviation in disor- 
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Figure 2. The top panel shows the structure factor (solid 
line: left y-axis) and stiffness (dashed line: right y-axis) vs. 
hopping for (jl/U = 0.7, zV/U = 1, pU = 20, and A/U = 
0.1 computed using MFT. The bottom panel shows the same 
as the top panel but computed with QMC on two different 
system sizes, L — 8 and 10. Both system sizes yield nearly 
indistinguishable results. The squares (circles) show results 
for the structure factor (stiffness) and the lines are guides 
to the eye. The inset uses finite size scaling at t/U = 0.05 
to show that L = 8 approximates the thermodynamic limit. 
Both the top and bottom panels show that both methods 
support overlapping order parameters (i.e., a supersolid) for 
weak disorder. 



der averaging. 

The bottom panel of Fig. [2] demonstrates the stability 
of the solid and the supersolid in the presence of quantum 
fluctuations and spatial disorder. Two different system 
sizes are plotted but yield indistinguishable results. We 
find that the structure factor and the stiffness have con- 
verged to the same value for L > 6 (see the inset). The 
identification of a supersolid in QMC (bottom panel of 
Fig. [2]) therefore qualitatively agrees with the same iden- 
tification in MFT (top panel) for weak disorder. 

Our parameter choices so far do not offer strong evi- 
dence for an intermediate disordered phase between the 
solid and the supersolid. Starting from the incompress- 
ible solid phase the theorem of inclusions [T3] implies a 
Griffiths-type scenario as we drive the system to a com- 
pressible state in infinite system sizes and at zero tem- 
perature. For the parameters considered here we expect 
that the disordered solid phase plays a role analogous to 
the Bose-glass in the disordered BH model. Thus the dis- 
ordered solid should always, in principle, appear between 
the solid and the supersolid (dashed line in Fig. [I]). But 
this region remains very difficult to probe because large 
V should considerably narrow the accessible parameter 
space for the disordered state. A similar narrowing is ob- 
served for the large U disordered BH model [H]. More 
importantly, our MFT ignores quantum fluctuations and 
the QMC is done on finite sized systems. We therefore 
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Figure 3. (Color online) The top panel shows the structure 
factor (solid line: left y-axis) and stiffness (dashed line: right 
y-axis) vs. disorder strength for fj,/U = 0.4, zV/U = 1, /3U = 
20, and t/U = 0.05 computed using MFT. The bottom panel 
shows the same but computed using QMC for L — 6 (left 
and right triangles), L = 8 (diamonds and up triangles), and 
L = 10 (circles and squares) where the lines are guides to 
the eye. Both methods imply that increasing just disorder 
strength can drive the solid into a supersolid. The insets show 
distribution histograms (Pd) for A/U = 0.6 and L — 10. 



expect very little evidence of the Griffiths-type region in 
our study. Our methods (QMC and MFT) do accurately 
probe regions far from this Griffiths-type region. We con- 
clude that the supersolid is stable because it appears far 
from the transition. 

We increase disorder to test the extent to which the 
solid (and supersolid) withstand spatial disorder in the 
chemical potential. To estimate the effects of increasing 
disorder on the solid phase we first increase A in our 
MFT. The top panel of Fig. [3] tracks the structure factor 
and the stiffness as we increase disorder starting from 
the solid. Here we see that disorder starts to trigger a 
supersolid near A/U ~ 0.2. The theorem of inclusions 
requires that the solid becomes a disordered phase for 
A = -E|/2 followed by another transition at A c > £|/2. 
Our MFT estimate for the onset of a supersolid at A/U ~ 
0.2 is consistent with this requirement because it is larger 
than a MFT estimate for the disappearance of the solid 
for A at E s g /2U « 0.15 with t/U = 0.05. 

We construct a mean field site-percolation picture to 
understand the disorder induced supersolid on a SC lat- 
tice. In the parameter regime considered here, the per- 
colating supersolid is a percolating superfluid coexisting 
with a solid. Consider a point in the solid phase on the t 
vs. /i phase diagram near the solid-supersolid transition 
in the A = limit. MFT treats each site individually. 
Increasing disorder will tend to move individual sites into 
a state nearby in phase diagram, i.e., the supersolid, if 
possible. For a uniform probability distribution of chem- 




Figure 4. (Color online) QMC results for a single disorder 
configuration for the parameters of Fig. [3] with L = 10 and 
A/U = 0.6 but for /3U = 50. We find solid order among all 
sites but we only plot a site i if it has large density fluctua- 
tions: (nf) — (rii)' 2 > 0.2. We draw all available bonds from i 
within the cube. The supersolid is allowed to percolate from 
one edge to another only if local density fluctuations per- 
mit local superfluid order (the off-diagonal matrix elements 
in Eq. (5f) along enough connected bonds. 



ical potentials in the interval \p — A, /i + A] a single site 
will have an energy in a specified range [fi + ei, /i + £2] 
with a probability P(ei < e < £2) = (ea — £i)/2A. Using 
the mcanficld boundaries of the solid phase and the site 
percolation threshold for a SC lattice (p c = 0.31 28 ), 
we estimate a critical disorder strength for percolation 
of the supersolid through the solid, ~ 0.13. This esti- 
mate is consistent with QMC calculations presented be- 
low. Percolation thus provides an intuitive picture for the 
onset of a supersolid with increasing disorder in SC lat- 
tices. We note that the non-zero stiffness and structure 
factor completely distinguish the percolating supersolid 
from the Bose glass and a percolating superfluid. 

We use QMC to include quantum fluctuations as we 
increase disorder starting from the solid phase. Figures [3] 
(bottom panel) and [4] show QMC results suggesting that 
quantum fluctuations allow disorder to trigger a percolat- 
ing supersolid from a solid. Large disorder strengths limit 
QMC calculations. For example, convergence for the 
data point at A/U = 0.8 and L = 10 in Fig. [3] required 
~ 10 4 CPU hours on R410 PowerEdge Dell servers. 

The top and bottom panels of Fig. [3] are qualitatively 
consistent. In the bottom panel we see that as we in- 
crease A: i) The solid remains robust for A/U < 0.2. 
ii) The stiffness has a slight, L dependent upturn for 
0.2 < A/U < 0.4. Here we expect a Griffiths-type phase, 
the disordered solid, in the thermodynamic limit and 
P — > 00. iii) The supersolid begins to percolate through 
the lattice for 0.4 < A/U < 0.6 where L independent 
data support nonzero values for both S v and ps- iv) A 
weak supersolid for A/U > 0.6 (where the percolating su- 
persolid order parameters may have L dependence) and 
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eventually a Bose-glass at large A. 

We have established a platform for the study of disor- 
dered supersolids. Our MFT and QMC results both show 
that supersolids are stable against weak spatial disorder 
on the SC lattice. The supcrsolid remains stable far from 
expected Griffiths-type regions in parameter space. We 
also find a striking transformation of the solid into a su- 
pcrsolid with increasing disorder strength. Further work 
will explore critical properties and map out the full QMC 
phase diagram. 
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